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Abstract Discrete-time random walks and their extensions are common tools for an¬ 
alyzing animal movement data. In these analyses, resolution of temporal discretiza¬ 
tion is a critical feature. Ideally, a model both mirrors the relevant temporal scale 
of the biological process of interest and matches the data sampling rate. Challenges 
arise when resolution of data is too coarse due to technological constraints, or when 
we wish to extrapolate results or compare results obtained from data with different 
resolutions. Drawing loosely on the concept of robustness in statistics, we propose a 
rigorous mathematical framework for studying movement models’ robustness against 
changes in temporal resolution. In this framework, we define varying levels of robust¬ 
ness as formal model properties, focusing on random walk models with spatially- 
explicit component. With the new framework, we can investigate whether models 
can validly be applied to data across varying temporal resolutions and how we can 
account for these different resolutions in statistical inference results. We apply the 
new framework to movement-based resource selection models, demonstrating both 
analytical and numerical calculations, as well as a Monte Carlo simulation approach. 
While exact robustness is rare, the concept of approximate robustness provides a 
promising new direction for analyzing movement models. 
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1 Introduction 


Major advances in tracking technology during the last decades have made large 
datasets of animal movement available to ecologists, and analyses of data have be¬ 
come widespread in ecology. These analyses have shed light on mechanisms that 
underly fundamental processes such as migration (Robinson et al||2009[[Costa et al 


|2012| ), navigation ( [Tsoar et al|201 lt[Benhamou et al|2Ql it , or home range behaviour 
and territoriality ( |Borger et al||2008( ~ Potts and Lewis ||2()T4[ |Giuggioli and Kenkre| 
2014| ). They have helped to identify conservation goals by revealing habitat prefer¬ 


ences and critical environmental features for populations ( [Sawyer et al|2009||Colcherc| 
et al|2010]|Ito et al|2Q13[|Masden et al|2Q12| ), as well as the role of important mutu- 


alistic interactions between mobile animals and immobile plants ( [Cortes and Uriarte 


2013[ Mueller et al|2Q14 1. These are only few of the many facets of movement ecol¬ 


ogy- 

Mathematical and statistical models provide a framework for studying movement 
( Schick et al|2008t[Smouse et al|201Q[[Langrock et al|2013[ ). When linking models 
to data, we can estimate model parameters and identify best-fitting models, thus in¬ 
ferring unknown quantities or mechanisms in movement behaviour. Although move¬ 
ment itself is a continuous process, many individual-based movement models treat 
time as a discrete variable, viewing movement as a series of locations in space, or 
equivalently as a series of steps ( [Turchin[[l998t [McClintock et al[[2014[ ). This may 
largely be ascribed to data being available in this format. Discrete-time models may 
thus be an intuitive first choice to describe a sampled movement path. However, there 
may be more reasons to use discrete-time models. The continuous movement path of 
an animal may consist of various behaviours at different scales ( [Johnson et al|2002l 
Benhamou[[2013[ ). Using a discrete-time model at the scale of interest allows us to 


focus on the behavioural mechanisms at that scale, while, for example, combining 
other unknown processes as stochastic effects. Also, the choice of time formulation 
in a movement model can have side effects that impact inference results. For exam¬ 
ple, McClintock et al ( 2014[ ) demonstrated that using a continuous-time Ornstein- 
Uhlenbeck process in a hierarchical model for identifying behavioural states led to 
difficulties discriminating between states, due to an inherent correlation between the 
variables step length and bearing in the Omstein-Uhlenbeck process. 

Many movement models are based on random walks ( [Turchin|1998[ Codling et al 


2008 [[McClintock et al 2014| ). From simple random walks that assume independently 


and identically distributed steps, we have moved to correlated random walks, which 
include directional persistence ( [Kareiva and Shigesada[[198^ , and biased random 
walks, which can, for example, be used to model centralizing tendencies or long-term 
directional goals ( [Benhamou|2006[[Borger et al|2008[ [McClintock et al|2012] ). Many 
animals live in heterogeneous environments, and the composition of the environ¬ 
ment and availability of resources influence movement decisions ( [Fortin et al|2005t 
[McPhee et al||2012l >. Therefore, another trend of random-walk extensions has left 
behind assumptions about environmental homogeneity in favour of spatially-explicit 
models that incorporate habitat effects on movement decisions ( [Rhodes et al[[2005[ 
Avgar et al|2013l|Potrs et al|2014| i. These models have an advantage over statistical 


resource-selection and step-selection functions ( Manly et al|2002[ Fortin et al|2005[ 
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Forester et al||2009] ) by allowing simultaneous estimation of movement parameters 


and environmental effects. 

When linking discrete-time models to data, the temporal resolution of the dis¬ 
cretization is a critical feature that must be chosen with care. Different time scales 
may come into play and need to be consolidated. On the one hand, a time scale is 
given by the biological process of interest. For example, we may be interested in in¬ 
ferring behavioural mechanisms of a movement process and thus need to consider 
the time scale at which these mechanisms are relevant. The discretization of a model 
should represent this scale appropriately. On the other hand, a different time scale 
may be given by the data collection rate. In practice, the sampling rate of data is 
subject to technological constraints. One of the major limitations of electronic tag¬ 
ging devices such as Argos or GPS tags is battery life, imposing a tradeoff between 
measurement rate and total deployment time ( |Ryan et n||2004[ Breed et n||2011 1. 
Also, to avoid a large noise to signal ratio, the time interval should be chosen so that 


measurement error relative to distance travelled during a time interval is small (Ryan 


et al|2004] ). For slow moving animals and depending on the accuracy of the tagging 


device, a minimum time interval of an hour may be necessary (Jerde and Visscher 


2005 [ ). Therefore, the resolution of the data may not always match the time scale of 


the behavioural process of interest. In this case, it becomes a challenge for a model 
to overcome the conflict. 

A related problem is that sampling rate can affect data analysis results ( |Codling| 
land Hill|[2005t |Rowcliffe et al|[MT^ |Postlethwaite and Dennis||20T^ . A common 
measure calculated from raw movement data is the total distance travelled, which 
can provide useful information about an animal’s energetic expenditures. It is well 
documented that this quantity is highly influenced by the sampling rate of the data 
( |Ryan et ^|2004t |Mills et ^|2006 jTanfema et al||2012l |Rowcliffe et al|[MT^ . A 
range of studies demonstrated that other fundamental movement characteristics vary 
with data sampling frequency as well, for example path sinuosity and apparent speed 
( [Codling and Hill|2005| ), movement rate and turning angle ( jPostlethwaite and Dennis j 
2013| ), and estimates of territory size ( Mills et ^[2006 ). One of the main problems 


underlying these effects is information loss when subsampling a movement path. This 
also impairs our capacity to correctly estimate behavioural states through hierarchical 
modelling approaches that have become widespread in movement analyses ( [Breed 


[et al[2Ql l][Rowcliffe et al 2012 ^ . These findings demonstrate that great care is needed 
when extrapolating movement analysis results beyond the temporal scale of a study. 
Comparisons of results may not be appropriate if the temporal resolution of the data 
varies too much, but it is unclear what constitutes ‘too much’. 

Despite the evidence of the extent of the problem, little is known about how 
to solve it. Previous approaches have been mainly empirical, using very fine scale 
data or synthetic data from simulations, which are subsampled at various resolutions. 
Movement characteristics calculated at these varying sampling rates are then com¬ 
pared to the values based on the full data, which represent the ‘true’ values. Some 
studies have fitted functions to the relationships of movement characteristics and sam¬ 
pling rate ( [Pepin et al|2004[ [Codling and Hill|2005t [Mills et al|2006[ ). These empir¬ 
ically obtained functions may be used to correct movement characteristics for sam¬ 
pling rate. While correction factors derived from movement data remain situation- 
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specific and cannot easily be applied across species ( Ryan et al|2004 Rowcliffe et al 


2012| ), we can obtain more general results by analyzing the effects of sampling rate 


at the level of the model ( [Codling and Hill|2(3Q5{[Rosser et al|2Q13] ). Often, important 
characteristics of movement can be well captured by models, and therefore analyzing 
the properties of models can provide more general insights. However, only few such 
studies exist. An approach to circumvent the problem of scale-dependent statistical 
inference has been taken by Fleming et al[ ( [2014| ), who use the semivariance function 
of a stochastic movement process to identify multiple movement modes acting at dif¬ 
ferent temporal scales. The method takes into account all possible time lags between 
observations. However, there are limitations as to the movement processes that can 
be included in this analysis ( [Fleming et al|2014| ). 


Here, we present a rigorous framework for studying how movement models react 
to changes in sampling rate, and we use this framework to analyze a class of mod¬ 
els based on random walks. With our analysis, we seek to understand whether, and 
how, models can help to compensate mismatching temporal scales between different 
data sets or between data and behavioural process of interest. Focusing on spatially- 
explicit random walks, we investigate whether there are models that can validly be 
applied to data with different temporal resolutions and how we can account for the 
differences in resolutions in our interpretation of statistical inference results. In par¬ 
ticular, we are interested in how model parameters, and their estimates, change as 
we decrease the temporal resolution. While estimates may change due to a shift in 
behavioural scale, which we always need to be aware of, we are interested in the 
changes that arise from the method, that is the model. Our framework is related to the 
statistical concept of robustness, which aims at safeguarding statistical procedures 
against violations of model assumptions ( [Hampel| 19861 [Huber and Ronchetti|2009[ ). 
Often, such violations refer to deviations from assumed probability distributions (e.g. 
Normal errors), which may result in outliers, misspecified relationships between re¬ 
sponse and explanatory variables in regression analyses, or violations of the common 
independence assumption. In this paper, we define robustness of movement models 
against changes in temporal discretization. In our framework, we treat robustness as 
a formal property of a model, namely the movement model. If a model has this prop¬ 
erty, it can be applied to data with varying resolutions. Additionally, while model 
parameters do not stay the same, they change systematically and can be translated 
between resolutions. 


Our paper is outlined as follows. In section we define what we mean by a 
movement model to be robust against changes in temporal resolution. We provide 
three different definitions, varying in their strength of conditions. In section we 
present different approaches how the definitions can be used to analyze robustness 
of movement models. Depending on models’ complexity and preexisting informa¬ 
tion, we can use formal analytical methods, numerical calculations, as well as Monte 
Carlo and simulation approaches. We use these approaches to examine robustness of 
spatially-explicit random walks and resource-selection models, and we summarize 
our findings in section In section we discuss the relevance of our robustness 
framework for statistical inference and also draw specific conclusions for spatially- 
explicit resource-selection models. 
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2 Robustness of Markovian movement models 

We consider movement models that are discrete-time Markov processes of the form 
(X^, f G r), where Xt G M? is an individual’s location and T = {0, t,2t, ...} is a set 
of regularly spaced times. This means that we assume that the time interval T > 0 
between two successive location measurements is fixed. Such data often arise from 
terrestrial animals fitted with GPS devices ( [Frair et al|2010| ). The time interval T of 
the model is usually specified by the resolution of the data. We denote the one-step 
transition density for the probability of moving from location y to x between times 
t — T and t by 0), where 0 G 0 is a vector of model parameters. This 

notation highlights that the transition density can be time-heterogeneous. 

We consider sub-models that consist of every ^th location of the original model 
for ^ G N. The transition density of the nth sub-model for the probability of mov¬ 
ing from location y to x between times t — nt and t is denoted by pt-nT,t {x\y,e)- 
compare Fig. A priori, the function can have an entirely different form 

than pt-T,t and may correspond to a different probability distribution. However, via 
the Chapman-Kolmogorov equation, the ^-step transition density can be written as a 
marginal density, 

Pt—nT^t{^^t\^t—nT') 

— • • • 5-^/—(^—fzT: T • • • 1)T5 (1) 

JR^x---xR^ 

where we marginalize over all intermediate locations visited between times t — nz 
and t. For simplicity, we use the general subscript ‘joint’ to denote any joint density 
of multiple locations. From the notation of the locations it is clear which joint density 
is meant. The Markov property of the model allows us to stepwise split up the joint 
density as follows 


Pjornti^t^^t—T’! • • • 

= Pt—T,t{^t\^t—T: ^)- ( 2 ) 

We can continue this until we obtain 

Pt—nT,t {^t \^t—nT 7 ^ ) 

r 1 

J X • • • X 

Therefore, we can use the one-step densities to calculate the ^-step density; compare 
Fig. For statistical inference, and thus for our robustness concept, the model pa¬ 
rameter vector 0 plays a crucial role. Although the ^-step density may belong to a 
different distribution than the one-step density, equation ^ justifies that we use the 
same parameter 0 in the notation of the ^-step density as in the one-step density. 

We define robustness in terms of the one-step and ^-step densities of a model. 
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Definition 1 (Robustness of degree n) Let ^ G N be finite. A movement model of 
the above type is robust of degree n if there exists an injective function gn.&^& 
such that 


Pt-n-zAAy^ = Pt-%,t{Ay^Sn{6)) for &\\teT andx,:y e (4) 

This definition requires that the ^-step densities are of the same functional form as the 
one-step transitions, where parameters of the model are appropriately transformed 
via the function gn. This means that if a model is robust, the ^th sub-model is in 
fact the same as the original model but with systematically adjusted parameters. The 
parameter transformation gn allows us to extrapolate the original parameter Q to the 
coarser temporal discretization of the ^the sub-model. Additionally, we can use the 
^th sub-model to infer the parameter 6 of the original model, because we can invert 
gn{6). Note, however, that this rests on the assumption that the original model defines 
the process of interest. If, instead we start at the coarser resolution, we would also 
need surjectivity of the function g^ to conclude the existence of the finer model. 

Robustness of degree n has important implications. Given a behavioural process 
of interest, described by a robust model with parameter 0, we can apply the model 
not only to data with matching temporal resolution T but also to coarser data with 
resolution nT (e.g. double time interval for n = 2). The parameter estimate Xj/ that 
we obtain from the coarser data is in fact an estimate of g„(0). From this, we can 
infer the value of 6 via 6 = 8nH¥)- Additionally, robustness allows us to compare 
studies pertaining to the same behavioural process but using data sets with different 
resolutions. If 0 is the estimate based on the finer data, it can be extrapolated to the 
coarser scale via the parameter transformation g„(0), for all degrees n for which the 
model is robust. 

Robustness as in Definition is a strong condition that we do not expect to hold 
but in few special cases of the density 6). However, equation may hold 

up to a function v{x,y), where v is a bounded function that could also depend on n 
or T. For practical applications, such approximate or asymptotic robustness may be 
sufficient. Therefore, we provide two additional definitions. 

Definition 2 (Asymptotic robustness of degree n) Let ^ G N be finite. A movement 
model of the above type is said to be asymptotically robust of degree n if there exists 
an injective function : 0 ^ 0 and a function v : x x M+ ^ M+ with the 

property v(jc,y; t) — 1 = ^(t) on x x IR+, such that 

Pt-nt,t{x\y^ = Pt-x,t{x\y^8n{Q)) v{x,y; t) for all r e T and x,y G (5) 

Here, ^ denotes the Landau symbol for the order of a function. If a model is asymp¬ 
totically robust, the ^-step densities are not exactly the same as the one-step densities, 
as was required in Definition However, the discrepancy between the densities is 
bounded by a function that is proportional to T. More precisely, for an asymptotically 
robust model we have 


1-Ct< 


Pt-n'[,t{x\y,6) 

Pt-T,t{^\y>Sn{B)) 


^ 1 4 ~ 


( 6 ) 
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for all X, y and 0, for some constant C > 0. Therefore, if the time interval T of the 
model is sufficiently small, the ^-step density will closely resemble the one-step den¬ 
sity with appropriately adjusted parameters. Asymptotic robustness of degree n im¬ 
plies that robustness of degree n is achieved as T ^ 0, that is when the time interval T 
approaches zero. 

In applications, the time interval T may not be chosen sufficiently small for Def¬ 
inition to be useful. Therefore, we give a variation of Definition in which the 
function v does not depend on T. 

Definition 3 (Approximate robustness of magnitude 5 and degree n) Let n eN 

be finite. A movement model of the above type is said to be approximately robust 
of magnitude 5 and degree n if there exists an injective function : 0 ^ 0 and a 
function v : x ^ IR+ with the property 0 < 1 — 5 < v(jc,y) < 1 + 5 for all x, 

y G for a 5 > 0, such that 

Pt-n-z,t{x\y^^) = Pt-%,M\y,gn{Q))v{x,y) iov aWt £T andj;,yeM^. (7) 


Analogously to equation condition 0 can be written as 

p,-n-z,t{x\y,e) 


1-5 < 


< 1 + 5. 


Pt-xMy^gn{B)) 

In fact, we may consider two different types of magnitudes. Setting 

Pt-nx,t{Ay,^) 


v{x,y) := 


Pt-x,t{x\y’Sn{B))' 


( 8 ) 


(9) 


this function depends a priori on the parameters, that is we have v(jr,y; 0 ), and the 
magnitude is 5^. If max^ 5^ exists, then this is the overall magnitude for the model 
with all possible parameter values. The magnitude determines how close ^-step densi¬ 
ties are to the parameter-adjusted one-step densities. If 5 is small, then the correction 
function v is close to one everywhere, and thus the ^-step density has similar values 
as the one-step density over its entire domain. 

Asymptotic and approximate robustness have similar implications for inference 
as robustness, but only approximately. The quality of the approximation depends on T 
or the magnitude 5. Suppose we wish to estimate parameters of a behavioural process 
that we formulate in a model. Suppose we consider the time interval T as suitable 
for the process. If the model is robust of degree n, we can use data not only at the 
matching scale but also at a coarser scale. For example, if the model is robust of 
degree 2, we can use data obtained at time interval 2 t. Because the model is also 
valid for the coarser scale, we can translate parameter estimates between the scales 
via the function If a model is asymptotically or approximately robust, the model 
is not exactly but still approximately valid for the coarser scale. To see this, consider 
the likelihood function 


Li{e\{X(i,Xr,X2x,-■■■,})= n Pt-x,t{Xt\Xt-x-,Q)- 


( 10 ) 
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If a model is robust of degree n, the likelihood for data at time interval nT is 

L„{e\{xa 

• • • 5 }) — 

t^{nT,(n+\)T,...} (11) 

• • • 5 })• 

If a model is asymptotically robust, we have instead 

U{gn{e)) . (1 -CT+ ^(T^)) < Ln{e) < U{gn{6)) • (1 + CT + (12) 

omitting the notation of the data, which is the same as in equation O- Analogously, 
for approximate robustness we have 

Li(g„(0)) • (1 -5 + ff{5^)) < Ln{e) < Li(g„(0)) • (1+C5 + (13) 

Therefore, if a model is asymptotically or approximately robust of degree n, we 
may loosely write Ln{0) ~ Li (g„(0)), that is the likelihood functions based on data 
at time interval T and on data at interval nT are approximately the same. Thus, if data 
at time interval T is not available, we can analyze data at time interval nT instead, 
using the likelihood Li of the original model. Parameter estimates obtained in this 
way can be translated to the scale T by using the inverse parameter transformation 
g~^. Such results from statistical inference based on Li may be close to results based 
on the correct L„, which may be difficult to compute. How close results are depends 
on the quality of the approximations in Definitions!^ andvia T or 5. For example, 
if a model is approximately robust with a very small magnitude 5, the likelihood Li 
will describe data at time interval nT almost as well as Ln. 


3 Analyzing spatially-explicit random walks 

We used the robustness definitions to analyze spatially-explicit random walk mod¬ 
els. These models merge general movement tendencies of an individual with deci¬ 
sions based on specific characteristics of locations, such as environmental features 
and available resources. We investigated how the models react when applied to data 
with increasingly coarser temporal resolution. 

Our robustness definitions have two key features. First, the one-step transition 
densities of the model and the ^-step densities of the sub-models need to have the 
same form. Second, model parameters, which are parameters of the densities, need 
to be transformed by a known function We can assume different approaches to 
investigate robustness properties of a model, depending on whether we have a can¬ 
didate for the parameter transformation gn or not. If prior knowledge allows us to 
investigate robustness for a given or hypothesized parameter transformation, we can 
calculate and compare the ^-step density Pt-nT.t W3’,0) and the parameter-adjusted 
one-step density pt-T:,t{Ay^Sn{^))- By showing equality of the two densities, we can 
verify robustness. For complex models, analytical calculations may be difficult, or 
even impossible. In these cases, we may resort to numerical calculations, especially 
when approximate robustness is sufficient. 
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In many situations, we may not know gn a priori, nor have any anticipation. Or, 
we may have tested robustness for a hypothesized parameter transformation but got 
poor results. In these cases, we need to establish some information on possible forms 
of the parameter transformation. Additionally, for complex models, numerical cal¬ 
culation of the high-dimensional integral required for the ^-step density (compare 
equation may become inaccurate. A solution is then to draw on the ideas of 
Monte Carlo sampling. Monte Carlo methods and simulations are useful when prob¬ 
ability densities are difficult to compute in closed-form but can conveniently be sam¬ 
pled from (e.g., [Robert and Casella 2000| ). In the following, we demonstrate both 
approaches for analyzing movement models’ robustness. 


3.1 Analytical and numerical approach 

Spatially-explicit random walks can be created by merging two elements in the tran¬ 
sition density of the model. One component is the general movement kernel (v;y), 
which can be the transition density of any standard random walk, describing the prob¬ 
ability that an individual takes a step from y to v if there were no environmental in¬ 
formation available. A second part of the model, given by the weighting function 
wq^{x), rates each possible step based on the location x. The transition densities of 
the full model takes the form 




k0^{x;y) we^jx) 

lRkei(z;y)w02(z)dz' 


(14) 


The integral in the denominator serves as a normalization constant. 

For simplicity, we restricted our analysis to the one-dimensional case, that is we 
assumed that Xt G M. We further focused on Gaussian kernels ke^{x;y) = ka{x\y), 
where ka{x\y) is a Gaussian density with mean y and standard deviation a. The 
weighting function wq^{x) was assumed to be positive everywhere to ensure that 
equation defines a density. In the following we simply use 0 for the parameter 
vector of the weighting function, or, when it is clear which parameters refer to the 
weighting function, we drop the subscript for the parameter in the notation of the 
weighting function entirely. 

Note that the transition density ( p^ does not depend on time explicitly. Still, as 
the individual moves through space over time, the centre location y of the kernel 
shifts. Although the kernel is a function of the distance ||x —y|| only, the weighting 
function adds a spatially explicit component. Therefore, unless the individual remains 
at the same location, the transition kernel effectively changes at every time step. In 
the following, we omit the time-related subscript in the notation of the density and 
simply write pi for the transition density and pn for the ^-step density. The time 
interval of the original process is always assumed to be T. The distinction between 
one-step and ^-step density is still important, because the ^-step density is in fact an 
integral over multiple one-step densities; compare equation 0. 

We investigated whether we could find weighting functions we (x) such that the 
model with transition density is robust, asymptotically robust or approximately 
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robust. We started by verifying Definition for simple cases of the weighting func¬ 
tion for a fixed parameter transformation As highlighted above, the parameter 
transformation is a key element, translating parameters between different temporal 
resolutions. For the parameter of the Gaussian movement kernel we obtained a 
candidate for the transformation based on the linearity of the Gaussian distribution. 
If we only consider the kernel ka, we have a simple random walk with normally dis¬ 
tributed steps between locations. The ^-step density ^ is then the density of a sum 
of n normally distributed random variables, which is again normal with standard de¬ 
viation ^/na. Therefore, we assumed that the transformation of the kernel’s standard 
deviation was given by gn{<y) = ^/na. For the parameters of the weighting function 
we assumed that they remain unaffected, that is gn{6) = 6. This is an ideal property 
for a weighting function, as it guarantees validity of inference results across different 
sampling rates without further translation. 

In a next step, we used the same parameter transformation g„(cr, 0) = (v^cr, 6 ) 
to establish conditions on the weighting function such that the model is asymptot¬ 
ically robust. For this, we assumed that the parameter of the kernel, the standard 
deviation, was infiuenced by the time interval T, that is cr = cr(T). This refiects that 
an individual may travel larger distances during longer time intervals. Because of 
the linearity of the Gaussian distribution, we assumed the relationship g{t) = \fxod, 
for some m > 0. For certain conditions on the weighting function, we verified Defini- 
tionj^analytically for the robustness degree ^ = 2 by calculating the function v(v,y; t) 
and placing bounds on it. 

As alternative to an analytical approach, we can calculate the ratio of two-step 
and one-step density numerically to see whether we can find a function v(v,y; t) ac¬ 
cording to Definitionfor the degree n = 2. Define S(t) := max;^:^; |v(v,y;T) — 1|. 
Note that since step densities depend on T through cr(T), we may equivalently con¬ 
sider 5(cr). If this is independent of the other parameters 0, we can obtain the bound 
on V as 5 := max^j 5((7), if this maximum exists. More generally, we can consider 
v(x,y, (7, 0) and calculate ^^(cr) := max^^j |k(x,y; cr, 0) — 1|. This ^^(cr) is the mag¬ 
nitude of approximate robustness (degree 2) for a model with a fixed weighting func¬ 
tion, including parameter values. An overall magnitude for the family of models con¬ 
sisting of the model for all parameter values can be obtained as 5 : = max^j ^ 

We demonstrate these two numerical approaches with an example weighting func¬ 
tion. 


3.2 Simulation approach 
3.2.7 Resource selection models 


Resource selection analyses link animal location data and environmental variables to 
understand animals’ space-use patterns in relation to their habitat. These studies pro¬ 
vide insight into species’ preferences or avoidance of habitat characteristics, which is 
important information for wildlife management and conservation purposes ( [Hebble- 
white and Merrill [200^[Latham et al|20TT|[Squires et al|20f3] ). Central methodologi¬ 


cal elements are resource selection functions (RSF) and resource selection probability 
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functions (RSPF), describing the probability of selection of certain units (e.g. pixels 
of land) by an organism based on environmental covariates ( [Manly et al|2002t Boyce 


own m a 


mere statistical framework ([Boyce et al|2002[ 

Courbin et al|2013[), incorporated into 

spatially-explicit models ([Rhodes et al[[2005[ 

Aarts et al|2011[), and become part of 


refer to [Lele et af| ( |2013| ) for details about the distinction of RSF and RSPF and use 
RSF as a general term for both concepts, unless otherwise stated. 

We include resource selection in the spatially-explicit random walk with transi¬ 
tion density 111 by letting the weighting function take the form of an RSF, wq (x) = 
W 0 (r(x)), where r(x) = (ri (x),..., r„(x)) is a vector of resource covariates at location 
X. Each rj is a function over space, representing resource covariates such as elevation, 
biomass measures, land cover type, and much more. The transition density becomes 


pi{x\y,a,d) 


k<,{x\y)w 0 {r{x)) 

kk(j{z-,y)we{r{z))dz' 


(15) 


In practice, geographical information is spatially discrete, and therefore the normal¬ 
izing integral in equation ( p3] ) becomes a sum over pixels, or cells, of land. Note that 
we still restrict our attention to one-dimensional models. 


The RSF can take various forms, and here we consider the two most commonly 
used ones ( [Manly et al|2002t Lele and Keim|2006] ), the exponential RSF, 


H'exp('‘W) =exp(/3-r(x)) 


(16) 


and the logistic function. 


wiog(r(x)) 


exp {a + P ■ r{x)) 
l+exp(a + /3-r(x))’ 


(17) 


The vector j3 comprises all selection parameters with respect to resource covariates 
r. A higher selection parameter means stronger selection with respect to the corre¬ 
sponding resource. In the logistic form, a is an intercept parameter, which can shift 
the inflection point of the logistic function away from zero. The inflection point is 
the point where the logistic function attains a value of 0.5, that is where the probabil¬ 
ity of selecting a resource is 50%. If the exponential form ( p^ is used, an intercept 
similarly to the one used in equation (TJ) is not identifiable, because it cancels in the 
definition of the transition density Therefore we have omitted it in equation ( p^ . 
The function wiog has range (0,1) and can therefore be used to describe probabilities. 
This means that this form can be used as RSPF, which for a given location y specifies 
the probability that an animal selects this location, given the covariate values of the 
location. In contrast, the exponential RSF can only specify values proportional to this 
probability, with unknown proportionality constant ( [Lele et al|2013[ ). 
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3.2.2 Sampling models and sub-models 


We examined the two models with weighting functions Wexp and wiog for their ro¬ 
bustness. Because the weighting functions depend on space through environmental 
information r they are highly non-linear, and therefore the transition densities are 
difficult to examine analytically. Sampling probability distributions is a convenient 
work around and has the additional advantage that we can control parameters and iso¬ 
late processes of interest. We thus simulated sample trajectories from the model with 
transition densities ( p~5] ). The joint density of a movement trajectory (vi,..., v//) G 
of length G N is given by 

N 

P]omt{xi,---,XN,B)= Pl{xi,6) Y[pi {x,\x,-i,6). (18) 

t=2 


Thus, we sampled successively from the transition densities to obtain a full movement 
trajectory. We obtained samples from the subprocess Xn = (xi ,...) consisting of 
every ^th location by subsampling the full trajectories. These subsamples represent 
samples from the model with transition densities being the ^-step densities |•, 0 ). 

Because the models rely on environmental data, we simulated resource land¬ 
scapes as realizations of Gaussian random fields with exponential covariance model 
( [Haran 2011[ Schlather et al 2013 1. This resulted in spatially correlated resource land¬ 
scapes, thus ensuring realism; compare Figurej^in Appendixj^ The sampled move¬ 
ment trajectories were based on these simulated landscapes. To avoid confounding 
effects and to keep results as clear as possible, we assumed that the weighting func¬ 
tion was based on only one resource r, thus we have WQ{r(x)). With the exponential 
covariance model, we assumed that the covariance of resource values at two different 
locations is given by 


Cov(r(x),r(j)) (19) 

where 5' affects the decrease of the spatial autocorrelation with increasing distance. 

We sampled trajectories for varying parameter values. We used a G {5,6,7} and 
P G {0.5,1,1.5,2} in all combinations. In the model with logistic RSF wiog, we fur¬ 
ther combined the values a G { —1,—0.5,0,0.5,1} with all other parameters. For 
each parameter combination, we sampled 16 trajectories for 15,000 time steps each; 
compare Fig. |9|10| in Appendix For each of the 16 trajectories, we used a differ¬ 
ent resource landscape, repeating the same set of resource landscapes across different 
parameter combinations. The 16 landscapes were generated with varying spatial auto¬ 
correlation, 5- ranging between 200-500. This led to a total of 192 sampled trajectories 
for the model with exponential RSF and 960 trajectories for the model with logistic 
RSF. We subsample every trajectory at levels ^ = 1,..., 15, leaving 1000 steps for the 
coarsest time series. The subsample for ^ = 1 is the original trajectory. 


3.2.3 Analyzing parameters 

While the simulated trajectories represent samples from the original model with tran¬ 
sition densities pi (• | •, 0 ), the subsamples of the full trajectories provide us with sam¬ 
ples from the sub-models with ^-step densities 0). To learn about the model’s 
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robustness properties, we need to test whether the subsamples reconcile with the 
parameter-adjusted one-step densities P\{'\'^gn{Q)) fof some parameter transforma¬ 
tion gn. For a given parameter transformation, we can achieve this by analyzing the 
fit of the model with transitions P\{'\'^gn{^)) with the subsamples. When gn is un¬ 
known, or when the fit for a hypothesized gn is poor, we first need to investigate the 
behaviour of the parameters under subsampling to see whether we can find a function 
gn as required by our robustness definitions. 

Here, we both tested a priori expectations on the parameter transformation and 
searched for better alternatives. We estimated parameters for all trajectories and their 
subsamples using maximum likelihood optimization. The likelihood function for the 
full trajectories is given in equation ( p^ . For subsamples, we applied the same model, 
although we did not know whether subsamples of trajectories followed the same 
(parameter-adjusted) process as full trajectories. We expected parameter estimates 
for the full trajectories to be close to the values that we used during the simulations. 
We call these the ‘true values’, although deviations in the simulations are possible, be¬ 
cause simulated trajectories are realizations of stochastic processes. Our main interest 
are parameter estimates for the subsamples. To distinguish estimates from underlying 
true parameters, we denote the estimate with a hat, e.g. a. Ideally, the parameters of 
the subsamples should follow some function g^((7, a, j3), and so should the estimates. 
To see whether such a function exists, we fitted non-linear regression models to the 
relationship of parameter estimates of subsamples and the subsampling amount n. 
For each parameter, we fitted two models. One model was more restrictive and repre¬ 
sented a priori expectations, whereas the other model had an additional free parameter 
that allowed more fiexibility for the parameter transformation. 

The general movement kernel k has one parameter, the standard deviation a of 
the Gaussian distribution. This kernel describes the general movement tendencies of 
the animal, and G influences the distance covered in each step. With increasing sub¬ 
sampling, the temporal resolution of the movement path becomes coarser, and we 
thus expected the standard deviation of the kernel to increase. Each step in a sub¬ 
sample is in fact the accumulated result of one or several steps in the full trajectory. 
If the kernel is the only force driving the movement, the linearity of the Gaussian 
distribution caused us to ex pect the standard deviation of the kernel to increase as 
^/nG\ compare section 3.1 With additional resource selection, however, there may 


be deviations from this behaviour. 

For the resource selection parameters a and j8, an ideal behaviour would be that 
they remain unaffected by the subsampling, analogously to our assumptions in sec¬ 
tion |3T] In our model, we assume that each step is influenced by the RSF. One of the 
underlying assumptions of a traditional RSF is that it gives weights to locations in¬ 
dependently of the values of other locations, which means each location is weighted 
by its present resource only, without consideration of alternative locations. There¬ 
fore, resource selection parameters should be independent of the temporal resolution 
of the data. However, within the spatially-explicit movement framework, resource 
selection always occurs in the context of the current location and the available sur¬ 
rounding area as defined by the general movement kernel. Therefore, a change in the 
movement kernel due to increased subsampling may be accompanied by a change in 
resource selection parameters. 
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We fitted the non-linear regression models to the parameter estimates separately 
for each parameter combination. This means that in each regression, we fitted esti¬ 
mates of 16 trajectories and their subsamples. Because of our previous considerations 
about the kernel parameter a, we assumed a power relationship between the estimate 
a and the subsampling amount n, stratified by trajectories. We chose the stratification 
because trajectories were simulated on different landscapes. Also, for the resource 
selection parameters, especially when their true values were close to zero, estimates 
could vary between being positive and negative. In these cases, the stratification al¬ 
lowed for flexibility. The model for the estimate of the ^th subsample of trajectory i 
is 

+ £, 1 < < 15, 1 < / < 16, (20) 

where the error term £ is normally distributed with mean zero and positive standard 
deviation (^. The maximum likelihood estimate of h should ideally be close to 0.5, 
however as noted above, it may deviate from this value because of resource-selection 
mechanisms. To test whether b differs from 0.5, we used model selection via AIC 
between the model in equation ( [20| ) and the model in which we fixed Z? = 0.5. 

Model choice for the resource selection parameters was less clear. Visual inspec¬ 
tion of the estimates, preliminary fits with varying models and inspection of residuals 
suggested a power law for the parameter j3 as well. We thus fitted the following 
model, 

Pi,n = Ad • 1 < ^ < 15, 1 < / < 16. (21) 

We compared the fit of this model with the model in which we assumed that subsam¬ 
pling does not change the estimate by setting b = 0 . 

For the intercept parameter a in the logistic form of the resource selection func¬ 
tion, we chose a linear model, 

^i,n = ^i,l1)£^ l</t<15, 1 < / < 16. (22) 

Inspection of residuals suggested that in some cases the relationship between a and 
n was non-linear. However, a power-law model or other non-linear relationships were 
not consistently more suitable either. Therefore we remained with the simpler, the 
linear, model, noting that this is a mainly illustrative analysis. 

3.2.4 Calculating approximate robustness 

To accompany the simulation analysis, we examined approximate robustness proper¬ 
ties of the two models with exponential and logistic RSF. We focused on approx¬ 
imate robustness of degree 2, and we tested the ideal parameter transformations 
g2{<y,P) = {V2g,P) and g 2 {o,a,P) = (\/2(7,a,j3) forwexp and wiog, respectively. 
We numerically calculated a magnitude 5 = max;^: 3 ;(|v(v,y) — 1|) for every possible 
scenario that we used in the previous section. This means that we calculated a magni¬ 
tude for each combination of the parameters cr, j3, and a (in case of the logistic RSF) 
and for each of the 16 simulated resource landscapes. We may therefore think of 5 
as 5((7,a,j8,/), for 1 < i < 16; compare Fig.We examined whether magnitudes 
were influenced by parameter values and specific characteristics of the landscapes. 
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such as their spatial autocorrelation and their overall variation Var(r(x)) over the spa¬ 
tial domain. We further calculated an overall maximum max^j ^ 5(cr, We 

compared results between the model with exponential RSF, Wgxp, and logistic RSF, 

f^log- 


4 Results 

4.1 Analytical and numerical results 

We found few special cases of weighting functions wq that, together with the Gaus¬ 
sian kernel ka, resulted in a robust movement model according to Definition 

The simplest case was a constant weighting function. Such a weighting function 
reduces equation Gl to the case of a homogeneous environment, where only general 
movement tendencies play a role, but no environmental information. The model is 
then a simple random walk with normally distributed steps between locations. Be¬ 
cause of the linearity of the normal distribution, the model is robust of degree n for 
all ^ G N for the assumed parameter transformation g^(cr) = y^cr; compare also 
Theoremfor parameters a = b = 0. 

A natural next step was to consider a linear weighting function. However, a lin¬ 
ear weighting function violates the assumption of being strictly positive everywhere. 
If in equation ( p^ the current location y is the point at which w becomes zero, the 
normalization integral vanishes. Also, equation can become negative and thus 
cease to be a valid density function. Still, we could draw on the linearity of the ex¬ 
pectation of a random variable to look into this further. The normalization constant in 
the transition density ( p^ can be viewed as an expectation of the form E{w{Z)) for 
a normally distributed random variable Z with mean y. Therefore, if the function w is 
linear, the normalization constant reduces to w(y). Equation Gl then becomes 

pi{x\y,a,d) = (23) 

The right-hand side of the equation is positive whenever x and y are either both neg¬ 
ative or both positive. If movement only occurs in the domain where the weighting 
function is positive the model is robustness within this domain. The details of the 
proof can be found in Appendix [A[ 

Theorem 1 (Linear weighting function) Let w he a linear function w{x) = ax Ah, 
for ( 2 ,/? G M. Let ^ (ZM.be the interval where w > 0. Lor the restricted domain J’, 
the movement model with transition densities G1 is robust of degree nfor all G N. 
The parameter transformation is given by gn{o^a^b) = (y/na^a^b). 

We found another special case to be given by an exponential weighting function. 
Here, no restriction on the domain is necessary. Again, see Appendix [A| for details of 
the proof. 

Theorem 2 (Exponential weighting function) Let w be an exponential function of 
the form w(x) = Ce^^^^ for C^a^b G M. Then the movement model with transition 
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densities in is robust of degree n for all n with parameter transformation 
gn{o,C,a,b) = {y/na,C,a,b). 

The above two Theorems show that it is possible to verify exact robustness with 
the ideal parameter transformation gn{o^Q) = {\/nG, 6 ) for certain weighting func¬ 
tions. However, the cases are very restrictive, and robustness will fail for many other, 
and especially more complex, weighting functions. 

We could additionally establish asymptotic robustness for more general condi¬ 
tions on the weighting function. The main result is summarized in the following the¬ 
orem. For a detailed proof of the theorem, see Appendix]^ 

Theorem 3 (Asymptotic robustness of degree 2) Let wq be continuous and bounded 
away from zero. Let we further be twice differentiable with bounded second deriva¬ 
tive. Then the model with transition densities ( p^ is asymptotically robust of degree 
2 with parameter transformation g 2 (G, 0 ) = (v2cr, 6 ). 

Thus, if the weighting function is well-behaved according to the theorem, we can 
place a bound on the factor by which the one- and two-step density vary; compare 
equation El'- This bound is of order T, such that the discrepancy between one- and 
two-step density decreases with the time interval. 

Example 1 (Asymptotic robustness of degree As a simple example, consider the 
weighting function w{x) = sin(av) + j3 for a > 0 and j8 > 1. The choice of j3 guar¬ 
antees that the weighting function is positive everywhere. The function w is bounded 
between 0 < j3 — 1 < w(v) < j3 + 1 for all v G M, and its second derivative is bounded 
by |w"(v)| = a?. Therefore, Theoremholds. 

The proof of Theorem is constructive in the sense that it provides us with a 
constant C for equation ([^ in terms of the bounds on w and w". However, this con¬ 
stant may be rather large and does not necessarily provide the closest bound on the 
function v. Therefore, it can be informative to calculate approximate robustness nu¬ 
merically. 

Example 2 (Approximate robustness of degree 2) We continue the above example 
with weighting function w{x) = sin(av) + j8 for a > 0 and j3 > 1. We calculated the 
function v(x,y; cr, a, j3) from Definition numerically, using different values of a 
and j8 (Fig.|^). From this, we obtained S(^ p{G) (Fig.j^), which is the magnitude of 
approximate robustness (degree 2) for the model with specific weighting function (i.e. 
with specific parameters); compare Fig.[^ In each case, after reaching a maximum 
the function vanishes for increasing G. Therefore it appears that we can find 5(^ ^ := 
maX(j (cr). The wavelength of the sine curve, determined by a, and the intercept 
P have different effects on the function 5q^ p (cr). While a shifts the curve, j8 changes 
the height of the peak (Fig.lsb). Therefore, it appears that 5(^ ^ is independent of a 
and decreases for larger j3. For the weighting function to be positive, j3 needs to be 
larger than one. For j3 = 1, the function 5(^ ^ has a maximum at one. From these 
considerations, we can conclude that ^ ^ = 1. This is the overall magnitude 

of approximate robustness (degree 2) for the family of weighting functions w{x) = 
sin(av) + j8, a>0, j3>l; compare Fig. As a word of caution, we note that we 
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only calculated 5q^ p for a fixed number of parameter values and only within finite 
intervals for v and y, and therefore results may be limited to these ranges. 

In the region where 5(a) peaks, the approximation of the parameter-adjusted 
one-step density pi{x\y, \/2cr, a, j3) to the actual two-step density P2{x\y, cr, a, j3) is 
only rough. However, for larger values of a, and independent of a and j3, the func¬ 
tion 5(^ p{a) seems to vanish, which means that the approximation is good and the 
discrepancy between two- and one-step densities may be neglected. From Theorem|^ 
we would have been able to conclude that 5^^j 3 (cr(T)) is bounded by Ct, for a con¬ 
stant C > 0, for all a > 0 and j3 > 1. As we can see from the steep initial slope of 
5 (^1^(a), especially for higher values of a, the constant C would need to be rather 
large (Fig. [^). The calculations of approximate robustness could additionally show 
that the bound on v(v,y) is in fact much smaller. 


4.2 Simulation results 

4.2.1 Results for parameter estimates 

When analyzing parameter estimates from the simulated trajectories and their sub¬ 
samples, we found a difference in the behaviour of parameters between the exponen¬ 
tial and the logistic form of the RSF. Generally, subsampling had less effect on the 
value of parameter estimates using the logistic form, and the behaviour of estimates 
agreed closer with our expectations. 

For both RSF, estimates G showed a good fit with the power-law model. When 
we used the exponential RSF, the estimated power h ranged from 0.45 to 0.5 for 
varying parameter combinations, thus deviating from expected behaviour for some 
parameter combinations (Fig.|^). For small selection parameter j3, the estimate a 
showed the expected increase as a^/n. With increasingly strong selection, i.e. higher 
value of P, estimates G became smaller with increased subsampling relative to the 
ideal relationship. An increase in G did not infiuence the fit other than leading to 
a larger residual standard error (^, which is to be expected because of the overall 
larger values of the dependant variable. In contrast, when using the logistic RSF, the 
estimated power b differed only very slightly from 0.5 and in some cases, the simpler 
model with fixed b was preferred by model selection right away (Fig.|^). 

The behaviour of the resource-selection parameter j8 also differed between expo¬ 
nential and logistic RSF. For the exponential RSF, j3 showed a clear increase with 
increased subsampling, fitted well by our power-law model (Fig. [^). The power b 
remained similar (ranging 0.105-0.124) across parameter combinations, increasing 
slightly with larger G (Fig.[^). For the logistic RSF, estimates j8 generally remained 
closer to the original values for ^ = 1 (Fig.|^,d). In most cases, model selection via 
AIC preferred the power-law model to the ideal constant relationship, however, the 
estimated values of the power b are small, with 53 out of 60 values being below 0.1 
(total range 0-0.156, with one exceptional negative value Z? = —0.041). There was a 
tendency of b to be smaller and more concentrated under stronger selection (Fig. 15). 

Estimates of the intercept a in the logistic RSF showed a slight decline with 
increased subsampling in most cases (Fig.|^. This decreasing trend existed no matter 
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whether a was positive, negative, or zero. In general, slopes of the linear fit were 
all close to zero (ranging -0.047-0.058), and in a few cases the null model with b = 
0 was chosen. We found a trend in the realized intercept values in the simulated 
trajectories. With stronger effect of selection (larger j3), the intercept estimate a of 
original trajectories (^ = 1) was stronger concentrated around the true underlying 
value, which subsequently lead to a stronger concentration of estimates of subsamples 
(Fig.|6](. 

4.2.2 Results about approximate robustness 

When comparing magnitudes 5(cr, of approximate robustness (degree 2) be¬ 

tween the two models with exponential and logistic RSF, we found lower magnitudes 
for the model with logistic function wiog. Magnitudes for the model with exponential 
RSF ranged between 0.067 and 1.82, whereas those for the model with logistic RSF 
ranged between 0.02 and 1.19. The 5% quantile, the median and the 0.95% quantile 
were [0.092,0.34,0.97] (exponential RSF) and [0.046,0.21,0.64] (logistic RSF). 

We found that especially the selection parameter j8 had a strong infiuence on 
magnitudes, higher values of j8 leading to higher magnitudes (Fig.|^. For the model 
with exponential RSF, there was a tendency that weighting functions whose underly¬ 
ing landscapes had higher variation Var(r(x)) lead to smaller magnitudes (Fig. [^). 
However, we did not find an effect of the parameter 5' that was used in the simulations 
to infiuence the spatial autocorrelation of the landscapes. The model with logistic 
RSF did not show such an effect of landscape variation. The logistic model had the 
additional intercept parameter a. We found that higher values of a tended to result 
in lower magnitudes (Fig.[^). 


5 Discussion 

We have proposed a new rigorous framework for analyzing movement models’ ca¬ 
pacities to compensate for varying temporal discretization of data. Our framework 
comprises three definitions of varying strength for robustness of discrete-time move¬ 
ment models. Generally, if a model is robust, it can overcome problems of mismatch¬ 
ing temporal scales between different data sets or between data and biological ques¬ 
tions. Because our robustness is a very strong condition that holds only for very few 
and generally more simple models, we have introduced the additional concepts of 
asymptotic and, most importantly, approximate robustness. While for many move¬ 
ment models it is difficult, or even impossible, to examine the transition densities 
and their marginals analytically, approximate robustness properties of a model can be 
calculated numerically also for analytically intractable models. Therefore, we believe 
that especially approximate robustness will prove a useful new concept for movement 
analyses. 

We have formulated our robustness definitions in terms of the transition densi¬ 
ties of Markov models, because these models are often fitted to movement data with 
likelihood-based methods of statistical inference. For the considered models, we can 
obtain the likelihood function by multiplying the transition densities of subsequent 
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steps. If a model is robust, the transition densities keep their functional form across 
varying temporal scales, and parameters are transformed via a well-defined func¬ 
tion gn. The likelihood function therefore remains the same but will yield different 
parameter estimates. However, if the function parameter transformation is known, es¬ 
timates from one scale can be translated to estimates at another scales. If a model is 
only approximately robust, the likelihood function will not remain exactly but at least 
approximately the same under a change of scale. Depending on the magnitude of the 
approximate robustness, the approximation of the likelihood function may be suffi¬ 
ciently good to allow parameter estimates to be reasonably comparable for different 
scales, especially if the difference in scales is small. 

Our concept of robustness for discrete-time movement models is related to the 
formal concept of robustness in statistics. Generally speaking, robust methods in 
statistics acknowledge that models are approximations to reality and seek to protect 
outcomes of statistical procedures (e.g. hypothesis testing, estimation) against devi¬ 
ations from the underlying model assumptions. Classic examples are the arithmetic 
mean and median as estimates of a population mean: while the median is robust 
against outliers the mean is not (e.g. |Hampel||1986| ). Often, robustness is viewed in 
the context of deviations from assumed probability distributions (distributional ro¬ 
bustness; e.g. Huber and Ronchetti||20Q9| ). For example, data may be contaminated 
by few observations with heavier tailed distribution than the majority of the observa¬ 
tions. In regression analyses, robustness may also relate to the homoscedasticity as¬ 
sumption or the functional form of the response function (Wiens 2000} [Wilcox|20 1 2| ) . 
Additionally, robustness has been considered when the assumption of independence 
is violated and instead observations are correlated ( |Hampel||1986t [Wiens and Zhou| 
1996| ). In our paper, we consider robustness in the context of discrete-time move¬ 


ment models with respect to assumptions about the temporal discretization. In view 
of statistical robustness, we study violations against the assumption that the temporal 
resolution of our movement model, a stochastic process, matches the resolution of 
the data, when in fact the data is only a subsample of the assumed process. 

There is also a difference between our robustness of movement models and the 
well-established robustness in statistics. In our framework, robustness is a direct prop¬ 
erty of a model. In contrast, classical robustness in statistics is defined for objects such 
as estimators, test-statistics, or more generally, functionals (real-valued functions of 
distributions) ( |Hampel|1971 1986| ). For the type of models we have considered here, 
parameter estimates cannot be obtained analytically but through numerical optimiza¬ 
tion of the likelihood function. The likelihood function is build by the model’s transi¬ 
tion densities, and thus we have defined robustness at a very basic level. A possibility 
for future research is to investigate whether some of the formal concepts of statistical 
robustness can be applied to our framework to add further insight. With our paper, 
we provide a new perspective for studying effects of temporal discretization of move¬ 
ment processes, and we hope to encourage further research. 

Our analytical investigations indicate that robustness is a rare property among 
movement models, especially when behavioural mechanisms such as resource selec¬ 
tion are added. Therefore, if we apply models to data without considering this issue, 
we are in danger of misinterpreting results and drawing erroneous conclusions. How¬ 
ever, our analysis also shows positive prospects with respect to approximate robust- 
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ness. Theorem suggests that in slowly varying environments that produce locally 
linear weighting functions we may find some robustness. Theorem and the fol¬ 
lowing examples show that certain smoothness and boundedness conditions on the 
weighting function can lead to approximate robustness. In addition, Example fur¬ 
ther demonstrates that approximate robustness can be investigated numerically on a 
case-by-case basis. We have illustrated this with a smooth weighting function w{x) 
that is a direct function of space. In data applications, an animal’s preferences for 
locations usually do not depend on space per se but rather through the type of habi¬ 
tat and available resources, and the weighting function will be less regular. In our 
simulation study, we have therefore presented a case with a more realistic model. 

While it is difficult to analyze the transition densities and resulting ^-step densi¬ 
ties with analytical calculations, we have demonstrated with the simulation approach 
how we can still investigate robustness properties of complex models. Sampling from 
probability distributions instead of calculating them is the key idea of Monte Carlo 
methods. We have used this method to examine sub-models that have the ^-step den¬ 
sities as transition densities. With this we obtained the parameter transformation gn. 
Our approach differs from previous studies that have used subsamples of fine-scale 
data to establish an empirical relationship between sampling interval and movement 
characteristics ( [Pepin et al|20Q4t [Ryan et al|2004[|Rowcliffe et al|2Q12[ ). When using 
data, it can be difficult to tease apart effects that result from the methodology and ef¬ 
fects of actual behavioural changes at different scales. Analyzing model properties as 
we have proposed here allows us to identify those effects of temporal discretization 
that are attributable to the methodology. 

In our demonstration of the simulation approach, we analyzed spatially-explicit 
resource selection models. These models have an advantage over traditional resource- 
selection and step-selection functions. In the traditional, regression-type approach, 
observed movement steps are compared to potential steps that are obtained separately 
from an empirical movement kernel ( [Fortin et~al|[2005[ [Forester et al[[2009[ ). In this 
approach, movement and resource-selection are treated independently, although it is 
highly likely that both infiuence each other. In contrast, when fitting the full random 
walk with resource selection to data by using the likelihood function ( p^ , we can 
simultaneously estimate parameters both of the general movement kernel and the 
weighting function, that is the RSF. 

In our analysis of the resource-selection model, we observed systematic trends 
in values of parameter estimates with changing temporal discretization of movement 
trajectories. The main purpose was not to analyze these relationships in full detail 
but to illustrate that they occur and thus must not be neglected. Comparing the expo¬ 
nential and logistic form of the spatially-explicit resource selection model, we found 
that estimates varied more with increased subsampling when the exponential RSF 
was used, compared to the logistic RSF. Using the exponential RSF, estimates of the 
kernel standard deviation a decreased with increased subsampling compared to the 
ideal relationship ^/na. On the other hand, using the logistic RSF, a followed the 
ideal relationship that would occur for a purely Gaussian process very closely, even 
under additional infiuence of resource selection. The estimated strength of resource 
selection, indicated by j8, increased with the subsampling amount. While this effect 
was strongly pronounced for the model with exponential RSF, it was only weak for 
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the logistic RSF. Therefore, if using the logistic RSF, one may expect to obtain similar 
inference results across varying temporal discretization. 

When we calculated the magnitudes of approximate robustness for the models 
used in the simulations, we found that those were in line with the results for the pa¬ 
rameter estimates. Overall, the model with logistic RSF had better robustness prop¬ 
erties than the model with exponential RSF. We also found a matching trend for the 
movement parameter cr with varying true values of j8. Estimates of cr were closer 
to the expected behaviour for weaker resource-selection parameters. This was also 
reflected in magnitudes of approximate robustness. If selection was weaker in the 
original model, the model exhibited better robustness properties. These results sug¬ 
gest that numerical calculations of approximate robustness can assist our expecta¬ 
tions about changes in parameter estimates. On the other hand, although parameter 
estimates of the weighting function showed a clear difference in behaviour when 
comparing between the exponential and logistic RSF, differences within one model 
between different parameter combinations were less clear. More analyses would be 
required to entangle more detailed effects. 

Overall, the results from the simulations suggest that depending on the resolution 
of movement data, we may misinterpret animals’ movement tendencies and also may 
overestimate resource selection effects. It is therefore important that we are aware 
of the changes to statistical inference that can arise merely from the methodology. 
Here, we have seen that changes in inference results were stronger for the resource 
selection model with exponential RSF compared to the logistic RSF. A possible ex¬ 
planation may be the additional intercept in the logistic RSF. With increased sub¬ 
sampling, estimates of a tended to decrease, possibly counteracting the increase in 
estimates j3. This could have led to more stability for the parameter cr of the general 
movement kernel. However, this may not explain why resource selection parameters 
generally varied less themselves compared to the exponential RSF. Another possi¬ 
bility is that the different form of the RSFs causes their different behaviour. While 
the exponential form of the RSF greatly enhances differences in landscape values, 
the logistic RSF is restricted to values in the interval (0,1). Theorem [^suggests that 
variation in the rate of change of the weighting function influences robustness prop¬ 
erties. Thus the logistic RSF may produce more stable inference results for varying 
temporal resolutions. Lele and Keim (2006) suggested several alternatives to the ex¬ 
ponential RSF. Our study case showed that the choice of resource selection functions 
can have implications for statistical inference and we encourage to choose resource 
selection functions more deliberately. 

With our study we have illustrated that the concept of the parameter transforma¬ 
tion gn can help to bridge the gap between different temporal resolutions of data. 
In the model with exponential RSF, we found that with increased subsampling esti¬ 
mates of the resource selection parameter j3 deviated strongly from the original val¬ 
ues. However, the increase in j3 could be fitted with a power-relationship. Thus, using 
the idea of Monte Carlo sampling, we were able to obtain a parameter transformation 
gn. Using such transformations when comparing results obtained from data with dif¬ 
ferent temporal resolutions could greatly improve our statistical inference, leading to 
a better understanding of movement behaviour. 
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Pt —2T,t —T (^t — T l^t — 2t 5 — — 



Pt-2rA^t\Xt-2T,0) 


Fig. 1 The second sub-model consists of every second location. The transition densities of the sub-model, 
which we refer to as 2-step densities, are the marginals over the two intermediate one-step densities of the 
original model 



Fig. 2 Steps for calculating the magnitude of approximate robustness of degree 2 for a given model, where 
G is the parameter of the movement kernel, and a and j8 are parameters of the weighting function. The 
one-step density pi can, for example, be equation d with the weighting function from Example]^ or the 
resource selection model d with weighting function d or d- When the resource selection model is 
used, the flowchart shows the calculation of the magnitude for one speciflc resource landscape r(x). When 
calculating an overall magnitude, practically we do this for a subset of the parameter space 
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Fig. 3 Panel a): Numerical calculation of the function v(x,j), which is the ratio of two-step den¬ 
sity O', a,j 8 ) and one-step density p^_T:^^(x|y,g 2 (o',a,j 8 )), for the weighting function w(x) = 

j 8 + sin(ax). Parameter values are cr = 1, a = 1, j 8 = 2. The function v(x,y) varies roughly between 0.72 
and 1.31. Panel b): Numerical calculation of 5(o') := maxxj |v(x,y;o') — 1| for the weighting function 
w(x) = j 8 + sin(ax) for varying values of a and j 8 . The parameter a, which determines the wavelength 
of the sine, shifts the curve 5(o') and varies the skewing, while retaining the height of the maximum. The 
parameter j 8 in contrast changes height of the maximum, which is the magnitude 5 of the approximate 
robustness 
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□ 


estimates 


Fitted model 

— | b=0.5 (fixed) 
b estimated 


Fig. 4 Values of cr against increasing subsampling amount n. Estimates d (gray dots) were fitted with a 
power-relationship, stratified by trajectories, and separately for several combinations of true parameter val¬ 
ues (cr, j8, and a for the model with logistic RSF). The power b was either fixed at 0.5 (ideal relationship; 
upper orange lines) or flexible and estimated (lower blue lines). The noted range of b refers to variation for 
different parameter combinations. Estimates and predictions are standardized by the corresponding true 
value. Panel a): Model with exponential RSF. With increasing value of j8, estimates a tended to increase 
less with subsampling compared to the ideal relationship. Panel b): Model with logistic RSF. The fitted 
power-relationship was very close to the ideal relationship, such that lines indicating the ideal relationship 
are overlaid by lines showing the fitted relationship 
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Fig. 5 Simulation results for the resource selection parameter j8 for the model with exponential RSF (pan¬ 
els a,b) and logistic RSF (panels c,d)- Panels a) and c): Estimates j8 (gray dots) for increasing subsampling 
amount n were fitted with a power-relationship, stratified by trajectories, and separately for several com¬ 
binations of true parameter values (cr, j8, and a for the model with logistic RSF). The power b was either 
fixed at zero, representing the assumption that resource-selection parameters do not change with changing 
temporal resolution (ideal relationship; straight orange lines), or flexible and estimated (curved blue lines). 
Estimates and predictions are standardized by the corresponding true value. In panel c), only estimates and 
predictions for a = 0, j8 = 1 are shown. Panel b): For the exponential RSF, the estimated power b was 
always above 0.1 and tended to increase with a. Panel d): For the logistic RSF, the estimated power b was 
mainly below 0.1 and tended to decrease and concentrate more for increasing j8 
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□ estimates 

Fitted 
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b estimated 


Fig. 6 For the model with logistic RSF, values of a against increasing subsampling amount n. Estimates 
were fitted with a linear relationship, stratified by trajectories, and separately for several combinations of 
true parameter values (cr, j8, and a for the model with logistic RSF). The slope b was either fixed at zero, 
representing the assumption that resource-selection parameters do not change with changing temporal 
resolution (ideal relationship; straight orange lines), or flexible and estimated (blue lines). Estimates and 
predictions are standardized by the corresponding true value and only shown for a = 0.5. The noted range 
of b refers to variation for different parameter combinations 
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Fig. 7 Magnitudes of approximate robustness for the study case models with resource selection. The plots 
depict 5 for varying values of a and selection parameter j8 (dots). Lines join values for the same landscape 
i, 1 < i < 16. Panel a): Magnitudes for the model exponential RSF. Values of 5 tend to be lower for 
landscapes with less variation Var(r(x)). Panel b): Magnitudes for the model with logistic RSF. Values of 
5 tend to be lower for higher values of the additional intercept parameter a 
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A Proofs of results about exact robustness 


Proof (Theorem^ First, note that for any standard deviation of the kernel, C7, the integral ka (y;x)w(y) dy 
reduces to the weighting function evaluated at the kernel’s mean, 

/ ka{y,x)w{y)dy = / ka{y,x){ay+ b)&y = / ka{y,x){a{y - x + x)+b)&y 
7m 7m 7m 

= {ax + b) / ka{y,x)dy + a / ka{y,x){y — x)diy = ax + b = w{x)^ (24) 

7m 7m 

because ka(-\y) is a Gaussian density integrating to 1 and with vanishing first central moment. If we 
consider w as a linear transformation of a Normally distributed random variable with mean x, then equa¬ 
tion (24} reflects a special case of Jensen’s inequality, in which equality holds. 

We now show robustness of degree n with parameter transformation gn(o,a,b) — {^a,a,b) by 
induction. For = 1, we have the trivial transformation g\{a^a^b) = {a^a^b), and there is nothing to 
show for robustness of degree 1. 

We assume that robustness or degree n holds, that is we have the relationship 

Pn{xn\xo,o,a,b) = pi{xn\xQ, y/uG,a,b). (25) 

for all x„,xo G M. For zi + 1, we use the Chapman-Kolmogorov equation and Markov property and obtain 

n n-\-\ 

Pn+\{xn^\\xQ,G,a,b)^ / \\ p\{xk\xk-\,G ,a,b) doc\.. .doCn 

= / pi(xn+\\xn,G,a,b)[ / TTpi(x^tfe-i,CT,fl,Z?)dxi...dx^-i dx„ 

7m \ 7m«-i j 


-- / pi{Xn+\\Xn,G ,a,b) Pn{Xn\xQ,G ,a,b)d 0 Cn 

7m 

: / p\(Xn^\\Xn,G ,a,b) p\{Xn\xQ,^G ,a,b)^n 

7m 


(26) 
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where the last step follows by induction. We can now insert the model’s step probabilities and use equa¬ 
tion (24} to further calculate, 


. I U\ f kaiXn+l-,Xn)w(Xn+l) k^^{Xn,Xo)w(Xn) 

Pn+i{xn+i\xo,a,a,b)= - \ f -?- n / n , 

Jr J^k^{y;xn)w{y)dy kk^^{y;xo)w{y)dy 
ka{Xn+i;Xn)w{Xn+l) k^^{Xn,XQ)w{Xn) 


-I 

Jr 


w{Xn) 




-dx„ 


VV’fa+l) 

w(xo) 


/ ka{Xn+\\Xn)k /^^{Xn\XQ)diZ. 

Jr ^ 


(27) 


Note that we have assumed that all movement steps are within the domain where the weighting function 
is positive. Since ka{xn^\\Xn) = ka{xn+\ — x„;0), the integral in the last expression is the convolution of 
two Gaussian densities with variances and no^ and with means 0 and xq, respectively. Because of 
the linearity of Gaussian random variables, this is again a Gaussian density with mean xq and variance 
(n+ l)c7^. Because equation ( |^ holds for the kernel with any standard deviation, we can rewrite the 
denominator as w(xo) = /]^k^/^^^(y;xo)w(y)<iy. Thus, 


Pn ^\{ xn ^\\ xQ , a , a , b ) 


k^a{Xn+\\Xo)w{Xn+l) 

lR^vk+TcT(y’^oMy)^y 


= /7i(Xn+i|xoO 


(28) 


□ 


Proof (Theorem^ We proceed analogously to the previous proof. The integral of weighting function and 
kernel with arbitrary standard deviation cr and mean x is here given by 


f ka(r,x)w(y)dy= f ka(r,x)Ce‘‘y+‘’dy 
JR JR 

C f ( {y-xf 


JTmo . 


+ ay + b]&y. 


By completing the square and using substitution u = (y — x — ao^) we obtain 


/ ka{y\x)w{y)&y-- 
Jr 


y/lno 




y—x—aa 

Vlo 


dy 


a^^ax+b 




/ exp (—V^crdu. 
7 m 


The final integral reduces to \/^(7, and therefore. 


/ 

7m 


ka{y\x)w{y)&y = Ce 


(29) 
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Again, we prove robustness of degree n by induction, using parameter transformation gn{oa^b) = 
{^G,C,a,b). In the induction step, we obtain, with help of equation (2^ , 

I I h\- [ ka{Xn^\\Xn)Ce^^^+^^^ k^^{Xn,XQ)Ce^^^^^ 

Pn+i Xn+iXQ,o,a, j^ka{y\Xn)Ce^y^^dy S^k^^{y\XQ)Ce^y+^dy 

_ r ka{xn+i;xn)Ce^^-+^+^ k^^(x^;xo) ^ 

Jr 


Ce- 

g^n+l 




-dXn 


/ /:a(x„+i;x„)k^^(x„;xo)dz 

JR ^ 


£^n+\ 

^ ^-2- +axQ 


fa +1 ;xo) Ce^^n+i +b 
lRk^lcy{r^^o)Ce^y+^dy 


= pi{Xn+l\xo,y 


(30) 

□ 


B Proof of result about asymptotic robustness 


To highlight the main steps necessary to prove Theorem]^ we establish a series of intermediate results. As 
a first step, we show that the 2-step transition density can be broken up into a product of the form ^ in 
Definition 12] 

Proposition 1 The 2-step transition density of model with transitions d can be written as 

P2{Xt\Xt-2T,(y,0) ^ P\{Xt\Xt-2x,V2o,e)-v{Xt,Xt-2x\^), (31) 


where the function v is given by 


v(Xt,Xt-2x-,T:) 


lRkV2aiy'^^)y^e{y)<^y 

lRk(j{y,x)w0{y)dy 




jz) 

lRkc >{ r , z ) w0 { y)dy 


(32) 


Note that v depends on t through cr. For later convenience, we define 

lRkf2(y(y'’^)^eiy)^y 


Q(x;t) := ■ 


lRka { r , x ) w0 { y)dy 


/(xi,X 2 ;t):= f kcU\\{xx+X 2 )) dz. 

Jr V2\ 2 J jRka{y;z)w0{y)dy 


(33) 

(34) 


Proof The proposition can be shown with a straightforward calculation. The 2-step transition density is 
given by 


Pl{Xt\Xt-2xi^,^) 


= / 

Jr 


ka{xt\z)W0{xt) ka{z\Xt-2x)y^e (^) 

K /M^g(}';^)w0(}')dy 4kcy(y;x^_2T)w0(y)dy 

yveM [. ( w / N 

' ka{xt;z)k(y{z;xt-2x) 


dz 


Ir ka {y,xt-2x)we (y) dy Jr 


f 

Jr 


W0{z) 


lRkair , z ) w0 ( y)dy 


dz. 


(35) 

(36) 

(37) 

















Robustness of discrete-time movement models 


33 


Tthe product of the two Gaussian densities in the integrand can be transformed as follows 
ka {xt\z)ka{z\Xt-2z) = _2t ) ^ ^ 

v2 

The two-step density therefore becomes 
P2Mxt-2T,(y,o) 

_ k^^{xt;xt-2r)w0{xt) r /I 

jRk(y{y;xt-2T)w0{y)dy Jr ^ V' 2 

The numerator of the first factor is the desired one-step density up to appropriate normalization. If we ex¬ 
tend by the required normalization constant k^^ (y',Xt- 2 T )^0 (>^) dy, we obtain equations and ( |^ . 

□ 

We are now left to show that the function v — 1 is in the order of T on its entire domain x M+. In 
particular, this means that for any fixed z*, the function v(xi,X 2 ;t*) — 1 is bounded on via cT* for a 
constant c. It turns out to be helpful to analyze v separately on M? x (0, Tq) and x [to,oo) for some Tq. 
Because the proof is simpler for large T, we present this result first. 

Lemma 1 Let w be continuous and bounded away from zero, that is there exist L and U such that 0<L< 
W < U for all X G M. Let w further be twice differentiable on M with |w''(x) | < Mfor some M and all 
X G M. For any Tq > 0, we have v(xi ,X 2 ,; t) — 1 = ^(t) on x [tq, o®). 

Proof Let Tq be a number away from zero and fixed. Our goal is to establish bounds on the functions Q and 
/, as defined in and and to use these to place a bound on v — 1. Because w is twice differentiable 
we can apply Taylor’s theorem to obtain a linear approximation for w using any point x G M, 

(}') = (x) + vr' (x) (y - x) + R{y ), (40) 

where R{y) is the remainder term. This leads to 

/ ka{y;x)w 0 {y)dy = W 0 {x) / ka{y;x)dy Fw^x) / (y-■^)dy + / ka{y,x)R{y)dy, (41) 

t/M tzM t/K. t/M 

where the first term on the RHS becomes we(x), because the kernel integrates to 1, and the integral in 
the second term is the first central moment of the kernel, hence vanishes. The remainder R{y), using the 
Lagrange form, is given by R{y) = ^ (y — x)^, for some ^ between X 2 andy. Since the second derivative 
of w is assumed to be globally bounded, we have \R{y) \ < y (y — x)^. We use this to place bounds on the 
third term, recognizing that the remaining integral /]^kcr(y;x) (y — x)^dy is the second central moment of 
the Gaussian kernel ka, which is given by its variance = co^z. Therefore, 

< [ ka{y;x)w 0 {y)dy <W 0 {x)-\-^co^z. (42) 

2 JR I 

In general, the lower bound can be arbitrarily close to zero, therefore we cannot simply invert this inequal¬ 
ity to obtain an estimate on the inverse of the integral. Instead, we use the bounds on w and again the fact 
Ir k( 7 (y',x) dy = I for any cr and any x G M to establish 

0<L< f ka{y;x)w 0 {y)dy <U, (43) 

Jr 

which can be inverted. Since inequalities and ID hold for any cr and any x G M, they allow us to 
place bounds on both Q and /. For Q, we obtain 

^ [w 0 {x)-Mco^z) < Q{x\z) < i (w 0 (x) FMco^z) (44) 

for all X G M, T G . We can avoid the dependency of the bounds on x by again invoking the bounds on 
w. 


(z;^(x^+x^_2t)). (38) 


^ [L-Mco^z) < Q(x) < ^ -FMco'^z). 


(45) 
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For the function /, we only make use of the bounds on w and inequality ( |^ and get 

L , , U 

0 < — </(xi,x 2 ;t) < — 


(46) 


for all XI ,X2 G M, T G We can now continue to calculate v — 1. An upper bound is immediately given 
by 

' L 2 

With only few more additional steps, we obtain a lower bound by simply drawing upon L<U, its squared 
version and its inverse, 


v(xi,x2;t) - 1 = 0(xi;t)/(xi,x2;t) - 1 < ^ 


(47) 


, , , . U^-L^ ML 2 U^-L^ MU 2 

■ Wxl^ 2 ;t) - 1 ) < 


Define C := ^r the Tq chosen up front. Then, 


(48) 


U^-L^ MU 2 


11 < —-^ ffT+ 

(49) 



- O Cto 

(50) 

f T \ 



(51) 


The product on the RHS is non-positive for t > Tq, and hence |v(xi ,X2; t) — 11 < Ct for all x [tq, o^). 


The bounds on Q and /, and thus v — 1, established in the preceding proof are not sufficient to conclude 
the result as t ^ 0, unless L = U, which is the trivial case of a constant weighting function. More suitable 
bounds, however, can be found if inequality can be inverted. This can be achieved by assuming t to 
be small enough. 

Lemma 2 Let w be continuous and bounded away from zero, that is there exist L and U such that 0 <L< 
vrg (x) < U for all x G M. Let w further be twice differentiable on M with |w''(x) | < Mfor some M and all 
X G M. Let To = T^. Then v(xi,X 2 ,;t) — 1 = ^(t) on x (0,To). 


Proof Here we develop bounds on Q and I such that both Q —I and / — 1 are in the order of T. Let T < Tq 
for To as defined in the lemma. Then the lower bound of equation \A2) is bounded away from zero, 

/ \ 47 2 / \ 47 2 / \ 47 2 7L , , 

T> we(x)- —CO To > W0(x)- —CO =we(x)-L> 0 . ( 52 ) 


Hence we can invert the inequality and obtain 
wq (x) —Mco^t 


W0 (x) + y 0)2^ 


< G(x;t) < 


Wq (x) +MCU^T 
W 0 (x) - y m^T 


(53) 


Note that the values in the numerators and denominators differ slightly because the variances of the kernel 
k in the numerator and denominator of Q differ by a factor of 2. 

Since 2wq (x) — Mco^z >2L — Mco^Tq > 0, we can conclude 


e(x;T)-l < 


Wq{x) +Mco^t — Wq{x) — y m^T 


Mco^t 


Mco^t 


WQ(x)-^;<jffT 2 wq{x)-M co^z ~ 2 L-Mco^To' 

for all X G M and t < To. Using 2wq (x) +Mft)^T > 2wq (x) > 2L, we similarly obtain, 

3M(0^z 


(54) 


-(G(x;t)-1 ) < 


/ 3M 2 

< ^^(0 T 


2wq (x) + Mco^ T 2L 


(55) 
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for all X G M and T < Tq. If we set Ci := max (^2 L-2a^TQ ’ ^ follows that \Q(x;z) — 1| < CiT on 

X (0,To). 

Using analogous arguments as before, we can fine an find an upper bound on I, 


/(xi,X 2 ;t)= [ ka(z;l-(xi-\-X 2 )) 
m V 2 \ 2 / 


W 0 (z) 


lRka{y;z)we(y)dy 


dz 


W 0 (z) 


- dz 


W0 (z)-f 

W0 (z)-f m^T+f 


dz 


W0 (z)-f co^T 

= 


fco^z 

2 dz 


< 1 + 

A lower bound is given by 




fm^T 
L- yCU^To 


dz = 1 + 


we(z)- f co^T 
Mco^t 

IL-Mco^Zo ’ 


/(xi,X 2 ;t)> f ko ( z ;^{ xi + X 2 )) - 2 

jm V 2 y ^^^(z) + f 


I - [ ka (z;^(xi+X 2 )) * 
JM V2 V 2 / - 


T-dz^l- 


Mco^z 

2L 


Setting C 2 := 2 l-mco'^tq ’ obtain |/(xi,X 2 ;t) — 1| < C 2 T on x (0 ,To). 

We can now estimate v — 1 as follows, 

|v(xi,X 2 ;T)-l| = |e,/,-l|<|0,-l||/,-l| + |0,-l| + |/,-l| 

< Cl C2T^ -f (Cl + C2) T < (Ci C2T0 + Cl + C2) T, 

for all xi ,X 2 G M and all t < Tq. 


(56) 

(57) 

(58) 

(59) 

(60) 

(61) 

(62) 


(63) 

(64) 

□ 


Lemmata[^and[^ together with proposition[^prove Theoremj^ 
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C Supplemental Figures 



Fig. 8 Four of the simulated resource landscapes used for sampling movement trajectories. The de¬ 
picted landscapes have been generated with spatial autocorrelation Cov(r(x),r(j)) = exp(^^) for 
= 200,300,400,500. We standardized landscapes to range within the interval (—3,3). At the boundaries, 
we set values to -3 to avoid movement close to the boundary and thus boundary effects in the transition 
densities due to the normalization constant. 
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Fig. 9 Four of the simulated trajectories from the model with exponential RSF. The trajectories were gen¬ 
erated using the parameter values (7 = 6 and j8 = 1. The underlying resource landscapes are the landscapes 
depicted in Fig.jS] in the same order. 
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X 


Fig. 10 Four of the simulated trajectories from the model with logistic RSF. The trajectories were gener¬ 
ated using the parameter values cr = 6 and j8 = 1 (same as in Fig.[^ and a = 0. The underlying resource 
landscapes are again the landscapes depicted in Fig.|^ in the same order. 



















